/*!  \file  GeSDD.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      06/28/2013 12:39:41 PM
 *   \copyright GNU Public License.
 */
#include "GeSDD.h"
#include "Matrix.h"
#include "setThreads.h"
extern "C" {
void dgesdd_(
const char   *  jobz, 
const int    *     m, // Rows of the input matrix A
const int    *     n, // Columns of the input matrix A
      double *     a, 
const int    *   lda, 
      double *     s,
      double *     u, 
const int    *   ldu, 
      double *    vt,
const int    *  ldvt, 
      double *  work, 
const int    * lwork, 
      int    * iwork, // interger array dimntions = 8 * min(M,N)  
      int    * info);
}

GeSDD::GeSDD(): FullSVD( "Lapack dgesdd" ) {
	aColMajor = true; aRowMajor = false;
	setCompileTime(__DATE__, __TIME__);	
	canBeParallel = true;
}
bool GeSDD::setNumberOfThreads(int nThreads_) {	
	assert(canBeParallel);
	SET_NUM_THREADS(nThreads_);
	nThreads = nThreads;
	return true;
}
void GeSDD::compute(Matrix & A, Matrix & U, Matrix & Sigma, Matrix & VT ) {
	reset();
	int m = A.rows(), n = A.cols(), lda = A.lda(), 
	ldu = U.lda(), ldvt = VT.lda();
	int info = 0, lwork = -1;
	double wkopt = 0;	double* work = 0;
	int iSize = 8 * (n < m ?  n : m);
	maxMemory = 0; extraMemory = 0;
	int * iwork = NULL; allocate(iwork, iSize);
	if (!iwork) return;
	char job = 'A';
	dgesdd_( &job,  &m, &n, A.getModData(), &lda,	Sigma.getModData(), 
				U.getModData(), &ldu, VT.getModData(), &ldvt, &wkopt, &lwork, 
				iwork, &info );
	lwork = (int) wkopt;	allocate(work, lwork);
	if (work) { /*  Compute SVD using D&Q*/
		dgesdd_( &job, &m, &n, A.getModData(), &lda,	Sigma.getModData(), 
					U.getModData(), &ldu, VT.getModData(), &ldvt, work, &lwork, 
					iwork, &info );
		delocate(work);	
	};
	delocate(iwork);	
};
REGISTER(GeSDD);
